Materials von Tag 1 und 2
Die Materials von Tag 1 & 2 könnt ihr hier herunterladen. Alle files sind auf unter https://bit.ly/merk039, die Dokumentation des heutigen Tag 3 findet ihr unter https://bit.ly/merk035.
Mein Plan für Tag III
- Wdh: Einfache Multi-Level Begriffe konzeptuell verstehen
- No Pooling, Complete Pooling, Shrinkage
- Intercept Only Model, Random Intercept Model, Random Slope Model, Random Intercept and Random Slope Model
- Intraklassenkorrelation
- Fixed Effects, Random Effect
- Unsere Beispieldaten
- Einfache Multi-Level Modelle parametrisieren und schätzen
Unsere Beispieldaten
Sleep Study
Im Paket {lme4} ist der
Datensatz
sleepstudy enthalten. Diesen Daten liegt ein Experiment
zugrunde indem dem Pbd Schlaf entzogen wurde Für uns wichtige Variablen
sind:
- Reaction: Durchschnittliche Reaktionszeit in einem Reaktionszeitexperiment in (ms)
- Days: Anzahl der Tage des Schlafentzugs (der Schlafentzug startete nach dem zweiten Tag)
- Subject: Personidentifier
library(lme4)
head(sleepstudy)
| Reaction | Days | Subject |
|---|---|---|
| 249.5600 | 0 | 308 |
| 258.7047 | 1 | 308 |
| 250.8006 | 2 | 308 |
| 321.4398 | 3 | 308 |
| 356.8519 | 4 | 308 |
| 414.6901 | 5 | 308 |
BilaKi Daten
data_bilaki_wide <-
read_sav("data/BilaKi_wide.sav")
head(data_bilaki_wide)
| ID | V_ID | Tgroup | m0_s1 | m1_s1 | m2_s1 |
|---|---|---|---|---|---|
| 1 | 1 | 1 | 2.166667 | 3.166667 | NA |
| 2 | 1 | 1 | 2.333333 | NA | 3.166667 |
| 3 | 1 | 1 | 3.000000 | 3.416667 | NA |
| 4 | 1 | 1 | 3.666667 | NA | NA |
| 5 | 1 | 1 | 3.250000 | 4.000000 | 4.500000 |
| 6 | 1 | 1 | 4.500000 | 3.416667 | 4.833333 |
Elke hat uns
Daten
aus ihrem Projekt mitgebracht. Diese liegen im wide-Format
vor und müssen zunächst einmal in das long-Format gebracht
werden. Das kann bequem mit dem im tidyverse enthaltenen
Paket {tidyr} geschehen:
data_bilaki_long <-
data_bilaki_wide %>%
pivot_longer(names_to = "variable",
values_to = "perf",
cols = c(m0_s1, m1_s1, m2_s1))
head(data_bilaki_long)
| ID | V_ID | Tgroup | variable | perf |
|---|---|---|---|---|
| 1 | 1 | 1 | m0_s1 | 2.166667 |
| 1 | 1 | 1 | m1_s1 | 3.166667 |
| 1 | 1 | 1 | m2_s1 | NA |
| 2 | 1 | 1 | m0_s1 | 2.333333 |
| 2 | 1 | 1 | m1_s1 | NA |
| 2 | 1 | 1 | m2_s1 | 3.166667 |
Dann muss man noch aus den Variablennamen die Info über den Messzeitpunkt entnehmen und eine nominale Zeitvariable erstellen.
data_bilaki_long <-
data_bilaki_long %>%
mutate(time = substr(variable, 2, 2),
time_factor = factor(time))
head(data_bilaki_long)
| ID | V_ID | Tgroup | variable | perf | time | time_factor |
|---|---|---|---|---|---|---|
| 1 | 1 | 1 | m0_s1 | 2.166667 | 0 | 0 |
| 1 | 1 | 1 | m1_s1 | 3.166667 | 1 | 1 |
| 1 | 1 | 1 | m2_s1 | NA | 2 | 2 |
| 2 | 1 | 1 | m0_s1 | 2.333333 | 0 | 0 |
| 2 | 1 | 1 | m1_s1 | NA | 1 | 1 |
| 2 | 1 | 1 | m2_s1 | 3.166667 | 2 | 2 |
Popularity Data (Hox et al., 2017)
Dieser Datensatz entstammt den einführenden Beispielen des gut lesbaren Lehrbuchs von Joop Hox et al (2017) Für uns wichtige Variablen sind:
- popular: Eine Likertskala zur Selbsteinschätzung der Beliebtheit
- extrav: Extraversion (Big Five)
- texp: Berufserfahrung der Lehrkraft
Import
data_popularity <- read_sav(file ="https://github.com/MultiLevelAnalysis/Datasets-third-edition-Multilevel-book/blob/master/chapter%202/popularity/SPSS/popular2.sav?raw=true")
head(data_popularity)
| pupil | class | extrav | sex | texp | popular | popteach | Zextrav | Zsex | Ztexp | Zpopular | Zpopteach | Cextrav | Ctexp | Csex |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 1 | 5 | 1 | 24 | 6.3 | 6 | -0.1703149 | 0.9888125 | 1.486153 | 0.8850133 | 0.6690561 | -0.215 | 9.737 | 0.5 |
| 2 | 1 | 7 | 0 | 24 | 4.9 | 5 | 1.4140098 | -1.0108084 | 1.486153 | -0.1276291 | -0.0430845 | 1.785 | 9.737 | -0.5 |
| 3 | 1 | 4 | 1 | 24 | 5.3 | 6 | -0.9624772 | 0.9888125 | 1.486153 | 0.1616973 | 0.6690561 | -1.215 | 9.737 | 0.5 |
| 4 | 1 | 3 | 1 | 24 | 4.7 | 5 | -1.7546396 | 0.9888125 | 1.486153 | -0.2722923 | -0.0430845 | -2.215 | 9.737 | 0.5 |
| 5 | 1 | 5 | 1 | 24 | 6.0 | 6 | -0.1703149 | 0.9888125 | 1.486153 | 0.6680185 | 0.6690561 | -0.215 | 9.737 | 0.5 |
| 6 | 1 | 4 | 0 | 24 | 4.7 | 5 | -0.9624772 | -1.0108084 | 1.486153 | -0.2722923 | -0.0430845 | -1.215 | 9.737 | -0.5 |
Wdh: Begriffe
AA: Elaboriert diese Begriffe in dem Ihr Sie auf eines der händisch skizzierten Datenbeispiele (siehe live-demo) anwendet.
Einfache Multi-Level Modelle parametrisieren und schätzen
Sleep Study
Die Mehrebenenstruktur verstehen
Fragestellung
Wie kann die Mehrebenenstruktur der
Sleep Study Datensatzesbeschrieben werden?
Lösungshinweis 1
Mit head(sleepstudy) erkennt man, dass der Datensatz nur
aus AV, UV, Subjectidentifier besteht. Also muss können maximal zwei
Ebenen im Datensatz vorkommen.
Bildet man ein Tabelle mit Days und Subject
kann man entscheiden ob eine perfect iode rimperfecte Hierarchie
vorliegt.
Lösungshinweis 2
Eine solche Tabelle erhält man mit
table(sleepstudy$Days, sleepstudy$Subject)
Lösung
Es handelt um eine perfekte Zweiebenenstruktur mit folgenden Nestgrößen und folgender Nestanzahl:
library(lme4) #enthält den Datensatz
# Nestanzahl
length(unique(sleepstudy$Subject))
## [1] 18
# Nestgrößen
sleepstudy %>%
group_by(Subject) %>%
summarize(Nestgröße = n())
| Subject | Nestgröße |
|---|---|
| 308 | 10 |
| 309 | 10 |
| 310 | 10 |
| 330 | 10 |
| 331 | 10 |
| 332 | 10 |
| 333 | 10 |
| 334 | 10 |
| 335 | 10 |
| 337 | 10 |
| 349 | 10 |
| 350 | 10 |
| 351 | 10 |
| 352 | 10 |
| 369 | 10 |
| 370 | 10 |
| 371 | 10 |
| 372 | 10 |
Intraklasenkorrelation schätzen und verstehen
Fragestellung
Der ICC hilft dabei empirisch zu erkennen wie stark sich die Mehrebenenstruktur in der abhängigen Variable ausdrückt. Dazu schätzt man ein Intercept-Only Modell und setzt die Intercept-Varianz ins Verhältnis zur Gesamtvarianz.
\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]}, \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \]
Lösungshinweis 1
Die Syntax für das Intercept-Only Modell lautet
lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)
Lösungshinweis 2
Richtig schnell geht es mit der Funktion icc() aus dem
Package {performance}.
Lösung
library(performance)
library(lme4)
mod01_sleepstudy <- lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)
icc(mod01_sleepstudy)
## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.395
## Conditional ICC: 0.395
Dieser ICC ist vglw. groß! Multi-Level Modelling ist also nicht nur theoretisch impliziert sonder auch empirisch.
Single Level Model schätzen
Fragestellung
Um später entscheiden zu können inwiefern Random Interept und Slope empirisch gerechtfertigt sind, kann man zunächst ein Single Level Model schätzen um dieses später mit den komplexeren Multi-Level Modellen zu schätzen.
Lösung
mod02_sleepstudy <- lm(Reaction ~ Days, data = sleepstudy)
summary(mod02_sleepstudy)
##
## Call:
## lm(formula = Reaction ~ Days, data = sleepstudy)
##
## Residuals:
## Min 1Q Median 3Q Max
## -110.848 -27.483 1.546 26.142 139.953
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 251.405 6.610 38.033 < 2e-16 ***
## Days 10.467 1.238 8.454 9.89e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 47.71 on 178 degrees of freedom
## Multiple R-squared: 0.2865, Adjusted R-squared: 0.2825
## F-statistic: 71.46 on 1 and 178 DF, p-value: 9.894e-15
Das entspricht einem Modell \[ \operatorname{Reaction} = {\color{#8cd000}{\beta_{0}}} + {\color{#8cd000}{\beta}}_{1}(\operatorname{Days}) + {\color{#8cd000}{\epsilon}} \] mit den Koeffizienten \[ \operatorname{\widehat{Reaction}} = 251.41 + 10.47(\operatorname{Days}) \]
Random Intercept Model schätzen und vergleichen
Fragestellung
Wie ändern sich die Koeffizienten des Modells, wenn man zulässt, dass Intercepts zwischen den Clustern variieren?
Lösungshinweis
Das Intercept kann ja als konstanter Prädiktor mit dem Wert 1
aufgefasst werden. Daher lautet die Syntax für das Random Intercept
(1|Subject) (lies: »das Intercept darf über die
Clustervariable Subject variieren«).
Lösung
mod03_sleepstudy <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
summary(mod03_sleepstudy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1786.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2257 -0.5529 0.0109 0.5188 4.2506
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 1378.2 37.12
## Residual 960.5 30.99
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.4051 9.7467 25.79
## Days 10.4673 0.8042 13.02
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.371
Das entspricht einem Modell
## Warning: Colorization of greek notation not currently implemented for merMod
## models
\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] mit den Koeffizienten
## Warning: Colorization of greek notation not currently implemented for merMod
## models
\[ \begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1}}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(0, 37.12 \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] Die Modelle vergleichen kann man mit der Syntax
tab_model(mod01_sleepstudy, mod02_sleepstudy, mod03_sleepstudy)
| Reaction | Reaction | Reaction | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 298.51 | 280.65 – 316.37 | <0.001 | 251.41 | 238.36 – 264.45 | <0.001 | 251.41 | 232.17 – 270.64 | <0.001 |
| Days | 10.47 | 8.02 – 12.91 | <0.001 | 10.47 | 8.88 – 12.05 | <0.001 | |||
| Random Effects | |||||||||
| σ2 | 1958.87 | 960.46 | |||||||
| τ00 | 1278.34 Subject | 1378.18 Subject | |||||||
| ICC | 0.39 | 0.59 | |||||||
| N | 18 Subject | 18 Subject | |||||||
| Observations | 180 | 180 | 180 | ||||||
| Marginal R2 / Conditional R2 | 0.000 / 0.395 | 0.286 / 0.282 | 0.280 / 0.704 | ||||||
Die Modelle gegeneinander Testen \(H_0:
\text{beide Modelle klären gleich viel Varianz auf}\) kann man
mit dem anova()-Befehl:
anova(mod03_sleepstudy, mod02_sleepstudy)
## refitting model(s) with ML (instead of REML)
| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| mod02_sleepstudy | 3 | 1906.293 | 1915.872 | -950.1465 | 1900.293 | NA | NA | NA |
| mod03_sleepstudy | 4 | 1802.079 | 1814.850 | -897.0393 | 1794.079 | 106.2144 | 1 | 0 |
Unterschied ist signifikant, also steigert das Random Intercept die Modellpassung überzufällig.
Random Intercept Random Slope Model schätzen und vergleichen
Fragestellung
Wie ändern sich die Koeffizienten des Modells, wenn man zulässt, dass Intercept und Slopes variieren?
Lösungshinweis
Die Syntax für ddie Random Effects lautet
(1 + Dayes|Subject) (lies: »das Intercept und die Slopes
dürfen über die Clustervariable Subject variieren«).
Lösung
mod04_sleepstudy <- lmer(Reaction ~ Days + (1 + Days|Subject), data = sleepstudy)
summary(mod03_sleepstudy)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
## Data: sleepstudy
##
## REML criterion at convergence: 1786.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2257 -0.5529 0.0109 0.5188 4.2506
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 1378.2 37.12
## Residual 960.5 30.99
## Number of obs: 180, groups: Subject, 18
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 251.4051 9.7467 25.79
## Days 10.4673 0.8042 13.02
##
## Correlation of Fixed Effects:
## (Intr)
## Days -0.371
Das entspricht einem Modell
## Warning: Colorization of greek notation not currently implemented for merMod
## models
\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1j[i]}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{j}} \\ &\mu_{\beta_{1j}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{1j}} \\ \rho_{\beta_{1j}\alpha_{j}} & \sigma^2_{\beta_{1j}} \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] mit den Koeffizienten
## Warning: Colorization of greek notation not currently implemented for merMod
## models
\[ \begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1j[i]}}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &0 \\ &0 \end{aligned} \end{array} \right) , \left( \begin{array}{cc} 24.74 & 0.07 \\ 0.07 & 5.92 \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] Die Modelle vergleicht man mit der Syntax
tab_model(mod01_sleepstudy, mod02_sleepstudy, mod03_sleepstudy, mod04_sleepstudy)
| Reaction | Reaction | Reaction | Reaction | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 298.51 | 280.65 – 316.37 | <0.001 | 251.41 | 238.36 – 264.45 | <0.001 | 251.41 | 232.17 – 270.64 | <0.001 | 251.41 | 237.94 – 264.87 | <0.001 |
| Days | 10.47 | 8.02 – 12.91 | <0.001 | 10.47 | 8.88 – 12.05 | <0.001 | 10.47 | 7.42 – 13.52 | <0.001 | |||
| Random Effects | ||||||||||||
| σ2 | 1958.87 | 960.46 | 654.94 | |||||||||
| τ00 | 1278.34 Subject | 1378.18 Subject | 612.10 Subject | |||||||||
| τ11 | 35.07 Subject.Days | |||||||||||
| ρ01 | 0.07 Subject | |||||||||||
| ICC | 0.39 | 0.59 | 0.72 | |||||||||
| N | 18 Subject | 18 Subject | 18 Subject | |||||||||
| Observations | 180 | 180 | 180 | 180 | ||||||||
| Marginal R2 / Conditional R2 | 0.000 / 0.395 | 0.286 / 0.282 | 0.280 / 0.704 | 0.279 / 0.799 | ||||||||
Die Modelle gegeneinander Testen \(H_0:
\text{beide Modelle klären gleich viel Varianz auf}\) kann man
mit dem anova()-Befehl:
anova(mod04_sleepstudy, mod03_sleepstudy)
## refitting model(s) with ML (instead of REML)
| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| mod03_sleepstudy | 4 | 1802.079 | 1814.850 | -897.0393 | 1794.079 | NA | NA | NA |
| mod04_sleepstudy | 6 | 1763.939 | 1783.097 | -875.9697 | 1751.939 | 42.1393 | 2 | 0 |
Unterschied ist signifikant, also steigert der zusätzliche Random Slope die Modellpassung überzufällig.
Multi-level Modelle in MPlus
## Version: 1.1.0
## We work hard to write this free software. Please help us get credit by citing:
##
## Hallquist, M. N. & Wiley, J. F. (2018). MplusAutomation: An R Package for Facilitating Large-Scale Latent Variable Analyses in Mplus. Structural Equation Modeling, 25, 621-638. doi: 10.1080/10705511.2017.1402334.
##
## -- see citation("MplusAutomation").
## When hashfilename = FALSE, writeData cannot be 'ifmissing', setting to 'always'
## The file(s)
## 'mod.dat'
## currently exist(s) and will be overwritten
MPlus ist statistisch ultra-mächtig aber im vgl. zu sehr umständlich, da es keine Data Wrangling Funktionen hat. Daher muss man seine Multi-Level Daten zunächst mit einer anderen Software aufbereiten und als .dat-File speichern, bevor man dieses wieder in MPlus einliest, in einem .inp-Textfile die Syntax speichert und ausführt, aufdass die Ergebnisse dann in einem wiederum gesondeten .out-Textfile gespeichert werden. Da es im Mplus keine Objektorientierung gibt, muss man dann Modellvergleiche (z.B. Devianztest) händisch anhand der beiden .out-Files machen.
`
Die Syntax sieht in MPlus wie folgt aus:
TITLE: Random Intercept Random Slope Modell (sleepstudy)
DATA: FILE = "MPlus/mod.dat";
VARIABLE: NAMES ARE reac days subj;
CLUSTER = subj;
WITHIN = days;
ANALYSIS: TYPE = TWOLEVEL RANDOM;
MODEL: %WITHIN%
s | reac ON days;
%BETWEEN%
reac s;
JASP
In JASP/jamovi bekommt ihr das Random Intercept Random Slope Modell,
wenn ihr den sleepstudy Datensatz als
.sav-File
einlest und dann die folgenden Einstellungen vornehmt:
Einstellungen für ein Random Intercept Random Slope Modell
oder gleich diese reproduzierbare JASP-Analyse öffnet.
Literatur
Hox, J. J., Moerbeek, M., & Schoot, R. van de. (2017). Multilevel analysis: Techniques and applications (Third edition). Routledge.
Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models (Bd. 1). Cambridge University Press.
Snijders, T. A., & Bosker, R. J. (2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling (2nd ed). Sage.